package KIMS;	

public class Jacobi{

	/**
	 * 雅可比算法
	 * 计算实对称矩阵的特征值、特征向量
	 * @param rc 相关系数实对称矩阵，计算后矩阵的对角线元素就是源矩阵的特征值
	 * @param originortho 首次作为一个单位矩阵，递归时作为当前正交矩阵相乘的结果
	 * @param ortho 
	 * @param eps 精度
	 * @return 计算成功返回true，失败返回false（若增加"迭代次数"条件，则可能出现返回false的情况，否则只返回true）
	 */
	public static boolean computeJacobi(double[][] rc,double[][] originortho,double[][] ortho,double eps){

		int num=rc.length;//方阵的阶
		double[][] matrix1;//源矩阵经过正交变换后的矩阵
		double temp=rc[0][0];//临时替换变量
		double[][] trortho;//正交矩阵的转置
		int rownum = 0;//保存非对角线元素最大值的行标
		int columnnum=0;//保存非对角线元素最大值的列标
		double sine;//正弦
		double cosine;//余弦
		double theta;

		//选取对称对阵中非对角线元素中绝对值最大者
		for(int i=0;i<num;i++)
			for(int j=0;j<num;j++)
				if(i!=j)
					if(rc[i][j]<temp){
						temp=rc[i][j];
						rownum=i;
						columnnum=j;
					}
		theta=Math.atan(2*temp/(rc[rownum][rownum]-rc[columnnum][columnnum]))/2;
		sine=Math.sin(theta);
		cosine=Math.cos(theta);

		//给正交矩阵赋值
		for(int i=0;i<num;i++)
			for(int j=0;j<num;j++){
				if(i==j){
					if(i==rownum||i==columnnum)
						ortho[i][j]=cosine;
					else
						ortho[i][j]=1;
				}
				else if(i==columnnum&&j==rownum)
					ortho[i][j]=sine;
				else if(i==rownum&&j==columnnum)
					ortho[i][j]=-sine;
				else
					ortho[i][j]=0;
			}
		trortho=transposeMatrix(ortho);//正交矩阵的转置矩阵
		matrix1=multiMatrix(trortho,rc,ortho);//经过正交变换
		ortho=multiMatrix(originortho,ortho);
		for(int i=0;i<num;i++)
			System.arraycopy(ortho[i], 0, originortho[i], 0, ortho[i].length);

		//计算精度
		while(!comparePrecision(matrix1,eps))
			computeJacobi(matrix1,originortho,ortho,eps);
		return true;
	}

	/**
	 * 多矩阵的乘法
	 * @param matrix 需要进行计算的二维矩阵的数组（按乘法计算公式从左向右依次排列）
	 * @return 乘法计算结果
	 */
	private static double[][] multiMatrix(double[][]... matrix){
		
		double[][] temp=null;
		
		if(matrix.length<=1){
			System.out.println("矩阵个数小于1，无法进行乘法运算");
			System.exit(-2);
		}
		else{
			temp=matrix[0];
			for(int i=1;i<matrix.length;i++){
				int n1r=temp.length;//第一个矩阵行数
				int n1c=temp[0].length;//第一个矩阵列数
				int n2r=matrix[i].length;//第二个矩阵行数
				int n2c=matrix[i][0].length;//第二个矩阵列数
				if(n1c!=n2r){
					System.out.println("第一个矩阵的列数与第二个矩阵的行数不相等，无法进行计算！");
					System.exit(-2);
				}
				double[][] result=new double[n1r][n2c];//存放计算结果的矩阵
				for(int m=0;m<n1r;m++)
					for(int n=0;n<n2c;n++)
						for(int k=0;k<n2r;k++)
							result[m][n]+=temp[m][k]*matrix[i][k][n];
				temp=result;
			}
		}
		return temp;
	}


	/**
	 * 求源矩阵的转置矩阵
	 * @param matrix 需要转置的源矩阵
	 * @return 源矩阵的转置矩阵
	 */
	private static double[][] transposeMatrix(double[][] matrix){

		int row=matrix.length;
		int column=matrix[0].length;
		double[][] transpose=new double[column][row];//初始化转置矩阵

		for(int i=0;i<row;i++)
			for(int j=0;j<column;j++)
				transpose[j][i]=matrix[i][j];
		return transpose;	
	}

	/**
	 * 比较所得结果是否符合精度
	 * @param matrix 需要比较精度的矩阵
	 * @param eps 给定精度
	 * @return 如果符合精度要求返回true，否则返回false
	 */
	private static boolean comparePrecision(double[][] matrix,double eps){

		int num=matrix.length;

		for(int i=0;i<num;i++)
			for(int j=0;j<num;j++)
				if(i!=j)
					if(matrix[i][j]>eps)
						return false;
		return true;
	}
}